\documentclass[a4paper,10pt]{article}
\usepackage[pdftex]{graphicx, subfigure}
\usepackage[T1]{fontenc}
\usepackage[slovene]{babel}
\usepackage[utf8x]{inputenc}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{hyperref}


\begin{document}
\newcommand{\prm}{r_{M}}
\newcommand{\prp}{r_{P}}
\newcommand{\pqm}{q_{M}}
\newcommand{\pqp}{q_{P}}
\newcommand{\ptt}{\tau}
\newcommand{\pmm}{m}
\newcommand{\pnn}{n}
\newcommand{\pkk}{k}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}

\title{A synthetic approach towards building a custom biological oscillator}
\date{\today}
\author{Martin Stražar}

\maketitle



\section*{Abstract}

In the following article, we present an alternative method towards predicting the behaviour of a biological oscillator. Using a system of two differential equations, we obtain the parameters for a custom oscillator (defined by oscillatory frequency and amplitude). That is achieved by the usage of an efficient genetic algorithm and numerical methods.
In \cite{scheper99}, a mathematical model for a single negative feedback loop was proposed. The parameter space was not explored in its entirety. In this work, we propose an alternative approach to explore the parameter space: given the desired behaviour of the oscillator, the parameters that produce it are calculated.

\section*{Introduction}

%TODO: citiraj
Circadian oscillators are a vital part of every living organism. Many life processes depend on the periodicity of the internal clocks, like sleep and wake up cycles [],  oxygen and food consumption [], etc.
Many different mathematical and logical models have been proposed \cite{lema01}. Using a model, composed of a delayed negative feedback loop, first proposed by Goodwin \cite{goodwin65}, we perform a parameter space search in order to determine the proper values to achieve the desired behaviour.

%TODO: omeni repressilatorje, preveč komplicirano
%TODO: dodaj metode za dosego končnega cilja realizacije
%TODO: poišči primere nadzora velikosti celice, interakcij med proteini in to. Naj se nanaša na konkreten članek z rezultati. 
%TODO: post-trankripcijska regulacija. SiRNA, nastanejo iz dsRNA z vnosom endonukleaze (Dicer). Potencialni nadzor nelinearnosti v sodelovanju.


\section*{Methods and modelling}

%TODO: Naštej, citiraj
\subsection*{Model definition}
So far, many mathematical models of biological oscillators have been proposed. A common property is the negative feedback loop in the synthesis of an observed protein. Without further assumptions, the negative loop itself leads into a stable steady state []. Assuming an additional delay between transcription and translation phases and nonlinearity in protein synthesis cascade, the system produces oscillations in accordance with parameter values. Due to their dynamic nature during the protein synthesis, many processes are difficult to quantify. Examples include protein phosphorylation [] or protein folding, which is its own field of research altogether []. Consequently, numerical approximations are used in order to get as close as possible to a formal description of observed behaviour.
For the purpose of our proof of concept, we decided to keep the model as simple as possible. Scheper et al. \cite{scheper99} propose a model of intracellular circadian oscillator, based on the following Hill equations:
      \beq \frac{dM}{dt} = \frac{r_{M}}{1+(P/k)^n}-q_{M} \cdot M(t) \label{dmdt} \eeq
      \beq \frac{dP}{dt} = r_{P} \cdot M(t-\tau)^m-q_{P} \cdot P(t) \label{dpdt} \eeq
where $M(t)$ and $P(t)$ represent the concentrations of mRNA and the resulting protein, repectively. The paramters in the two equations are described in Table \ref{tab:partab}. 

\begin{table}[ht!]
    \caption{Parameters of the system an their initial values.}
    \label{tab:partab}
    \begin{tabular}{|l|l|l|}
        \hline
        Parameter label & Meaning   & Initial value                                  \\ \hline
        $r_{M}$              & Rate of production (mRNA)                  & 1.0\\ 
        $r_{P}$              & Rate of production (protein)               & 1.0\\
        $q_{M}$              & Rate of degradation (mRNA)                 & 1.0\\ 
        $q_{P}$              & Rate of degradation (protein)              & 1.0\\ 
        $\tau$               & Transcription-translation delay            & 1.0\\ 
        $k$                  & Scaling constant                           & 1.0\\ 
        $m$                  & nonlinearity in protein synthesis cascade  & 3.0\\  
        $n$                  & Hill coeficient                            & 2.0\\
        \hline
    \end{tabular}
\end{table}

Depending on the parameter values, the system can produce oscillations or end up in a stable steady state. In \cite{scheper99}, a parameter subspace search was performed, where not all the parameters in the equation were considered. In this work, we look at the problem from another angle: given the oscillator behaviour, we derive the parameter values to produce it.
Since it would be impossible to find all the solutions of the two equations analitically, we make use of a customized genetic algorithm to derive the proper parameter values. The input to the algorithm are the arguments of a sine wave, which describes the desired behaviour: 
  \beq y_{wave}(t) = A_{wave} \cdot sin(\omega_{wave} t) + \frac{A_{wave}}{2}, \eeq
where $A_{wave}$ and $\omega_{wave}$ stand for amplitude and frequency, respectively. In order to make the problem solution as general as possible, all quantities in the system are to be interpreted on an arbitrary unit scale.

\subsection*{Finding the parameter values}
An oscillator object is treated as a member of a population (individual). It is uniquely defined by the values of the eight parameters in table \ref{tab:partab}. Given the time interval and time step, the values of the equations \ref{dmdt} and \ref{dpdt} are computed. From the obtained numerical data, the following properties are derived: a flag, whether the oscillator is oscillating or not ($B_{osc}$), its amplitude ($A_{osc}$) and frequency ($\omega_{osc}$), and the rise and fall times ($t_{osc/fall}$ and $t_{osc/rise}$). These are used as metrics to evaluate the fitness function, which is the core of the selection process. The metrics are derived as described in table \ref{tab:metrics}.
\begin{table}[ht!]
    \label{tab:metrics}
    \begin{tabular}{|l|l|l|}
        \hline
        Shorthand & Metric & Evaluation \\ \hline
        $M_{A}$ & $Amplitude$                   & $|A_{wave} - A_{osc}|$                               \\ 
        $M_{\omega}$ & $Frequency$              & $|\omega_{wave} - \omega_{osc}|$                     \\ 
        $M_{sym}$ & $Symmetry$                  & $|t_{osc/rise} - t_{osc/fall}|$                     \\ 
        \hline
    \end{tabular}
\end{table}

%TODO: dodaj footnote
%TODO: popravi enačbe
%TODO: dodaj if else za osciliranje ali ne
The fitness function is then defined as follows:
\beq F_{osc} = \frac{M_{A} + M_{\omega} + M_{sym}}{3}\eeq


The initial population of oscillators is derived by adding small random changes in the initial values of parameters for each oscillator. The fitness function is evaluated for each oscillator and the best (the one with the smallest values of the fitness funtion) are preserved in the next generation. Additionaly, new oscillators are derived as offspring from the existing generation. The closer an oscillator is to the desired behaviour, the more offspring it will produce (the probability to get closer to the desired behaviour is larger). The number of offspring for an individual oscillator is calculated as follows:
\beq F_{mean} = E(\sum_{osc} F_{osc}) \eeq
\beq P_{rep}  = \{osc;  F_{osc} <= F_{mean}\} \eeq
\beq O_{total} = 1 + |P_{rep}|  \eeq
\beq O_{osc} = \frac{F_{osc}}{\sum_{osc \epsilon P_{rep}}} \cdot O_{total} \label{oshare} \eeq

where $F_{mean}$ is the mean fitness of the population, $P_{rep}$ is the part of the population to be reproduced, $O_{total}$ is the total number of offspring produced for the current generation, of which each individual oscillator gets its share by calculating equation \ref{oshare}. The child oscillator is derived from parent with a pretedermined probability for mutation in any of the definitive eight parameter values.
 In such manner, the procedure is repeated from generation to generation for a predefined number of times. In each generation new individuals are derived and the best are selected. According to \cite{man99}, such system will monotonically converge towards global maxima and sucessfully avoid local minima, which was the case in our tests as well. A more detailed description of the algorithm is given in the Appendix.

\subsection*{Improving the solution}
During the development phase, the initial solution was optimized in the following way. The simulation of the life of a population is repeated for more times, where the best solution at the end of one generation is the input to produce the next. We repeat such step for three times, where in each step the fitness function is modified, such that only a subgroup of the metrics is evaluated. A single step is repeated until the desired (user defined) accuracy for the given subgroup of the metrics is achieved. Steps are performed as shown in table \ref{steptable}. 
\begin{table}[ht!]
    \label{steptable}
    \begin{tabular}{|l|l|l|}
        \hline
        Step & Evaluated metrics  \\ \hline
        $1$              & Amplitude                             \\
        $2$              & Frequency                              \\
        $3$              & Amplitude, Frequency, Symmetry          \\ 
        \hline
    \end{tabular}
\end{table}

This modification allows us to find the solution more quickly, escape local maxima or to bias one of the metrics.


\section*{Results}
%TODO: slike
%run6
The algorithm was tested for various input behaviours. To reproduce the behaviour of \cite{scheper99}, with input amplitude of $16\ nM$ and wavelength of $24\ h$, we got perfectly symetric stable oscillations as close as $15.72\ nM$ amplitude ($0.017\ \%$ deviation) and  $24.14\ h$ ($0.0062 \%$ deviation). The sampling time was $1000\ h$ and sampling step was $6\ min$. The resulting parameters are shown in table \ref{tab:partab2}, column Test 1, and differ from the ones proposed in \cite{scheper99}. To achieve the result, 90 generations of oscillators were analyzed with maximum population size of 3200. The frequency and amplitude convergence graphs are show in figures [] and [].

%run1
In accordance with the tendency to achieve faster frequencies and tunable amplitudes, we tested the system for other input behaviours. In the following example, the input amplitude was $8.0\ nM$ and wavelength of $12\ h$. The end system resulted in perfectly symetric stable oscillations of $7.18\ nM$ ($10.23\ \%$ deviation) and $12.31\ h$ wavelength ($0.025\ \%$ deviation) and the resulting parameters are shown in table \ref{tab:partab2}, column Test 2.

\begin{table}[ht!]
    \label{tab:partab2}
    \begin{tabular}{|l|l|l|l|}
        \hline
        Parameter label & Meaning   & Test $1$ & Test $2$                               \\ \hline
        $r_{M}$              & Rate of production (mRNA)                  & 1.976 & 2.986\\ 
        $r_{P}$              & Rate of production (protein)               & 1.666 & 2.731\\
        $q_{M}$              & Rate of degradation (mRNA)                 & 0.536 & 0.739\\ 
        $q_{P}$              & Rate of degradation (protein)              & 0.996 & 1.134\\
        $\tau$               & Transcription-translation delay            & 9.236 & 4.262\\
        $k$                  & Scaling constant                           & 2.857 & 0.801\\
        $m$                  & nonlinearity in protein synthesis cascade  & 6.147 & 6.428\\
        $n$                  & Hill coeficient                            & 2.211 & 3.856\\
        \hline
    \end{tabular}
\end{table}

\begin{figure}[ht!]
     \label{fig:run}
     \begin{center}
        \subfigure[]{%
            \label{fig:run6}
            \includegraphics[width=0.45\textwidth]{img/run6.png}
        }
        \subfigure[]{%
            \label{fig:run6cmp}
           \includegraphics[width=0.45\textwidth]{img/run1.png}
        }\\ %  ------- End of the first row ----------------------%
        \subfigure[]{%
            \label{fig:run1}
            \includegraphics[width=0.45\textwidth]{img/run6cmp.png}
        }
        \subfigure[]{%
            \label{ffig:run1cmp}
            \includegraphics[width=0.45\textwidth]{img/run1cmp.png}
        }
    \end{center}
    \caption{%
      \ref{fig:run6} Test 1: concentrations of protein (full line) and mRNA (dashed line). \ref{fig:run1} Test 2: concentrations of protein (full line) and mRNA (dashed line). \ref{fig:run6cmp} Test 1: Comparison with protein concentration (full line) and target behaviour sine wave (dashed line).
     }%
\end{figure}


\begin{figure}[ht!]
     \label{fig:conv}
     \begin{center}
        \subfigure[]{%
            \label{fig:amp6}
            \includegraphics[width=0.45\textwidth]{img/run6amp.png}
        }
        \subfigure[]{%
           \label{fig:amp1}
           \includegraphics[width=0.45\textwidth]{img/run6frq.png}
        }\\ %  ------- End of the first row ----------------------%
        \subfigure[]{%
            \label{fig:frq6}
            \includegraphics[width=0.45\textwidth]{img/run1amp.png}
        }
        \subfigure[]{%
            \label{fig:frq1}
            \includegraphics[width=0.45\textwidth]{img/run1frq.png}
        }
    \end{center}
    \caption{%
       \ref{fig:amp6} Test 1: Amplitude value convergence. 
       \ref{fig:frq6} Test 1: Frequency value convergence. 
       \ref{fig:amp1} Test 2: Amplitude value convergence.
       \ref{fig:frq1} Test 2: Frequency value convergence.  
     }%
\end{figure}

%TODO: citiraj
The results show that there are multiple points in parameter space that produce similar behaviours for a given system, as well as that with a given desired behaviour, the required parameter values can be derived. Parameter values can be controlled to some degree (using promoters with various strengths [], proteins and operator sites with different binding rates [], post-transcriptional regulation to modify the delay [], etc.), allowing us to engineer arbitrary structures with desired behaviours.




\section*{Conclusions}

The advances in synthetic biology, genetic engineering and bioinformatics will allow researchers to shift towards a synthetic approach in buildic biological logic components. The growing registries of DNA parts (promoters, regulator binding sites, coding sequences etc.) help us to predict the properties of the systems that are to be built (example: http://www.biobricks.org). 
A number of health conditions require treatment in a predetermined and widely distributed intervals of time. For instance, to facilitate problems caused by insomnia, melatonine is taken after in defined periods of time \cite{paul04}. The insulin levels oscillate in periods of 5-15 minutes, where irregularities in the insulin oscillations are the main cause of diabetes type 1 \cite{porksen02}. To treat these and similar conditions, systems with predetermined oscillation period are ought to be built. By knowing the behaviour of parts and using proper modelling and prediction algorithms, the development is significantly faster.


\section*{References}
\begin{thebibliography}{9}
\bibitem{scheper99}
   Scheper T., Klinkenberg D., Pennartz C. and van Pelt J.,
   \emph{A Mathematical Model for the Intracellular Circadian Rhythm Generator} 
   The Journal of Neuroscience,
   1 January 1999, 
   19(1):40-47

\bibitem{man99}
  K. F. Man, K. S. Tang and S. Kwong,
  \emph{Genetic algorithms: concepts and design},
  Springer, 
  1999, 
	1-85233-072-4

\bibitem{paul04}
  Paul M. A., Gray G., Sardana T. M., Pigeau R. A.,
  \emph{Melatonin and zopiclone as facilitators of early circadian sleep in operational air transport crews}
  Aviat Space Environ Med.,
  2004, 75(5):439-43.

\bibitem{porksen02}
    Pørksen N., Hollingdal M., Juhl C., Butler P.,Veldhuis J. D. and Schmitz O,
    \emph{Pulsatile Insulin Secretion: Detection, Regulation, and Role in Diabetes},
    Diabetes, 
    February 2002, 
    vol. 51 no. suppl:245-254

\bibitem{lema01}
Lema M., Echave J., and Golombek D.,
\emph{(Too Many) Mathematical Models of Circadian Clocks (?)}
Biological Rhythm Research
2001, 
Vol. 32, No. 2, pp. 285–298

\bibitem{goodwin65}
Goodwin, B. C. 
\emph{Oscillatory behavior in enzymatic control processes}
Advances inEnzyme Regulation, 1965, 3:425-428.


\end{thebibliography}

\end{document}

